##REPLICATION CODE for Hazlett and Parente (2022)##

# List of packages for session
.packages = c("ggplot2", "sp", "maptools","remotes","sensemakr")

# Install CRAN packages (if not already installed)
.inst = .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])

# Load packages into session 
lapply(.packages, require, character.only=TRUE)

#Load data from sensemakr package
dat = colombia

########################
##VIOLENCE (Main Text)##
########################

##Table 2 (main text)

# The super naive model 1:
violence1.out   = lm(yes_vote ~
                        fat_2011to2015_gtd, data = dat)

round(summary(violence1.out)$coef,2)
summary(violence1.out)$df

#Sensitivity (model 1)
sense.out1 = sensemakr(violence1.out, treatment = "fat_2011to2015_gtd", kd = 1)

signif(sense.out1$sensitivity_stats$r2yd.x,2)
signif(sense.out1$sensitivity_stats$rv_q,2)
signif(sense.out1$sensitivity_stats$rv_qa,2)


# Model 2:
violence2.out = lm(yes_vote ~ fat_2001to2005_gtd + fat_2006to2010_gtd + fat_2011to2015_gtd + total_eligible + santos10 + gdppc , data = dat)

round(summary(violence2.out)$coef,2)
summary(violence2.out)$df

# Sensitivity (model 2)
sense.out2 = sensemakr(violence2.out, treatment = "fat_2011to2015_gtd", benchmark="santos10", kd = 1)

signif(sense.out2$sensitivity_stats$r2yd.x,2)
signif(sense.out2$sensitivity_stats$rv_q,2)
round(sense.out2$sensitivity_stats$rv_qa,3)



##Figure 1 (main text) 
plot(sense.out2, lim = 0.5, label.bump.x = 0.03, label.bump.y = 0.01)

##Sensitivity from existing studies

#Tellez 2019
signif(sensemakr:::robustness_value.numeric(t_statistic = 3.14, dof = 4200),2)
signif(sensemakr:::robustness_value.numeric(t_statistic = 3.14, dof = 4200,alpha=0.05),2)

#Pechenkina and Gamboa 2019
signif(sensemakr:::robustness_value.numeric(t_statistic = 3.25, dof = 1008),2)
signif(sensemakr:::robustness_value.numeric(t_statistic = 3.25, dof = 1008,alpha=0.05),2)

#####################################
##POLITICAL AFFILIATION (Main Text)##
#####################################

##Table 3 (main text)
outcome.santos2  = lm(yes_vote ~ santos14 + fat_2010to2013 + elev + gdppc + pop13, data = dat)
signif(coef(summary(outcome.santos2))[,1],2)
round(coef(summary(outcome.santos2))[,2],2)
signif(coef(summary(outcome.santos2))[,3],3)
summary(outcome.santos2)$df

sense.santos.out = sensemakr(outcome.santos2, treatment = "santos14", benchmark=c("gdppc","elev"), kd = 3)
signif(sense.santos.out$sensitivity_stats$r2yd.x,2)
signif(sense.santos.out$sensitivity_stats$rv_q,2)
round(sense.santos.out$sensitivity_stats$rv_qa,3)

##sensitivity from existing studies

#Krause 2017
signif(sensemakr:::robustness_value.numeric(t_statistic = 45, dof = 1069),2)

##############
##APPENDIX D##
##############

##Figure A1 (appendix)

remotes::install_github("nebulae-co/colmaps")
library(colmaps)

mapdat = read.csv("mapdata.csv")

#adding extra 0 to municipal codes for two departments to ensure match between colmaps and data
mapdat$id2 = ifelse(mapdat$department=="Antioquia",sprintf("%05d",mapdat$id),
                 ifelse(mapdat$department=="Atlantico",sprintf("%05d",mapdat$id),mapdat$id))

#Violence panel (A)

#create sum of fatalities from 2011-2015
mapdat$gtdsum = mapdat$gtd2011f+mapdat$gtd2012f+mapdat$gtd2013f+mapdat$gtd2014f+mapdat$gtd2015f

#log to shrink distribution
mapdat$log_gtd = log(mapdat$gtdsum+.001)

#transformation to create more granularity on continuous scale
mapdat$gtd_def = sign(mapdat$log_gtd)*abs(mapdat$log_gtd)^(1/7)

#plot
gtdmap = colmap(municipios,data=mapdat,var="gtd_def",data_id="id2")+
  scale_fill_continuous(low="#fee0d2",high="#de2d26",name="Exposure to \n         violence",breaks=c(1,-1),
                        labels=c("High","Low"),guide=guide_colourbar(direction="horizontal",title.position="bottom"))+
  theme(legend.position=c(0.28,0.16))

plot(gtdmap)

#Santos 2014 vote panel (B)
santosmap = colmap(municipios,data=mapdat,var="santos14",data_id="id2")+
  scale_fill_continuous(low = "#deebf7", high = "#3182bd",name="Support for Santos",breaks=c(75,25),
                        labels=c("High","Low"),guide=guide_colourbar(direction="horizontal",title.position="bottom"))+
  theme(legend.position=c(0.28,0.17))

plot(santosmap)

#Referendum panel (C)
votemap = colmap(municipios,data=mapdat,var="yes_vote",data_id="id2")+
  scale_fill_continuous(low = "#efedf5", high = "#3f007d",name="Support for peace",breaks=c(100,0),
                        labels=c("High","Low"),guide=guide_colourbar(direction="horizontal",title.position="bottom"))+
  theme(legend.position=c(0.28,0.17))
plot(votemap)

##############
##APPENDIX E##
##############

##Sensitivity from existing studies (appendix)

#Tellez 2019
sensemakr:::robustness_value.numeric(t_statistic = 3.14, dof = 4200)
sensemakr:::robustness_value.numeric(t_statistic = 3.14, dof = 4200,alpha=0.05)
sensemakr:::partial_r2.numeric(t_statistic = 3.14, dof=4200)
sensitivity_stats(estimate=0.22, se = 0.07, dof=4200,q=0.25)


#Pechenkina and Gamboa 2019
sensemakr:::robustness_value.numeric(t_statistic = 3.25, dof = 1008)
sensemakr:::robustness_value.numeric(t_statistic = 3.25, dof = 1008,alpha=0.05)
sensemakr:::partial_r2.numeric(t_statistic = 3.25, dof=1008)

#Krause 2017
sensemakr:::robustness_value.numeric(t_statistic = 45, dof = 1069)
sensitivity_stats(estimate=0.622, se = 1.2*0.0138, dof=1069) #making SE 20% higher

##############
##APPENDIX F##
##############

##Figure A2 (appendix)
sense.santos.out = sensemakr(outcome.santos2, treatment = "santos14", benchmark=c("gdppc","elev"), kd = 3)


plot(sense.santos.out,label.bump.x = 0.05,label.bump.y = 0,list.par=list(mar=c(4,4,1,1),pty="m"),asp=1,lim=0.9,lim.y=0.4)

##Figure A3 (appendix)
sense.santos.out.nobm = sensemakr(outcome.santos2, treatment = "santos14")

plot(sense.santos.out.nobm, type="extreme", lim=.9, r2yz.dx=c(1,.5,.3),legend.bty="l")

##Model using Santos 2010 as treatment instead of Santos 2014
outcome.santos3  = lm(yes_vote ~ santos10 + fat_2010to2013 + elev + gdppc + pop13, data = dat)
sensitivity_stats(outcome.santos3,treatment="santos10")

